#include "auxiliar.h"

// internal product
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN:
// a - vector
// b - vertor
// OUT
// internal product
double dot(const double *a, const double *b)
{
	return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; 
}

// compute face's area
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN:
// face - elemete object
// OUT
// area
double F_area(pFace face)
{
	int i;
	vector<pVertex> vertex;
	double xyz[3][3], n[3];

	M_GetVertices(face,vertex);	
	for (i=0;i<3;i++) V_coord( vertex[i], xyz[i] );
	
	double a[] = { xyz[1][0]-xyz[0][0], xyz[1][1]-xyz[0][1], xyz[1][2]-xyz[0][2] } ;
	double b[] = { xyz[2][0]-xyz[0][0], xyz[2][1]-xyz[0][1], xyz[2][2]-xyz[0][2] } ;

	computeCrossProduct(a,b,n) ;
		
	return .5*norm(n);
}

double norm(const double *n)
{
	return sqrt( pow(n[0],2) + pow(n[1],2) + pow(n[2],2) );
}

// swap edge's nodes ids
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
void swap(int *IJ)
{
	int temp = IJ[0];
	IJ[0] = IJ[1];
	IJ[1] = temp;
}

// multiply a tensor K by a vector vec
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN:
// tensor - a matrix 2x2 (written as a vector: double tensor[4])
// vec - vector 2x1 
// OUT:
// mult = tensor * vec
void multTensorVector2D(const double *tensor, const double *vec, double *mult)
{
	mult[0] = tensor[0]*vec[0]+tensor[1]*vec[1];
	mult[1] = tensor[2]*vec[0]+tensor[3]*vec[1];	
}

void swap(double &a, double &b)
{
	double temp = a;
	a = b;
	b = temp;
}

// converts a string read from file to a double
// =========================================================================================
// IN:
// string
// OUT:
// string converted
double strToDouble(string &str)
{
	return strtod(getSubString(str), 0);
}

// converts a string read from file to an interger
// =========================================================================================
// IN:
// string
// OUT:
// string converted
int strToInteger(string &str)
{
	return atoi(getSubString(str));
}

// filters string the numeric part - internal use
// =========================================================================================
// IN:
// string
// OUT:
// char for the numeric part
const char* getSubString(string &str)
{
	string::size_type loc = str.find( "=", 0 );
	string numberstr = str.substr(loc+1, str.size()-loc);	
	return numberstr.c_str();
}

void getEdgesData(pEdge e, double *a)
{
	char *str[4] = {"a_00", "a_01", "a_10", "a_11"};
	for (int i=0; i<4; i++) EN_getDataDbl(e,MD_lookupMeshDataId(str[i]),&a[i]);
}

void swapVertex( vector<pVertex> &vertex)
{	
	pVertex tmpVertex = vertex[0];
	vertex[0] = vertex[1];
	vertex[1] = tmpVertex;
}

bool isEdgeOnBdry(pMesh theMesh, pEdge edge)
{
	return (E_numFaces(edge)==1 && M_numRemoteCopies(theMesh,edge)==0)?true:false;
}

double getExactSolution(pVertex vertex)
{
	return .0;
}

void getFCenter(pFace face, double *center)
{
	vector<pVertex> vertex;
	M_GetVertices(face,vertex);
	getFCenter(vertex[0],vertex[1],vertex[2],center);
	vertex.clear();
}

void getFCenter(pVertex v1, pVertex v2, pVertex v3, double *center)
{
	int i;
	double coords1[3], coords2[3], coords3[3];
	for (i=0; i<3; i++)
	{
		coords1[i] = coords2[i] = coords3[i] = .0;
	}
	V_coord(v1,coords1);
	V_coord(v2,coords2);
	V_coord(v3,coords3);

	for (i=0; i<3; i++)
		center[i] = (coords1[i] + coords2[i] + coords3[i])/3.0;
}

// Make a vector given to nodes
// =========================================================================================
// IN:
// A,B - pointers for node coordenates
// OUT:
// V - pointer for vector
void makeVector(const double *A, const double *B, double *v)
{
	for (int i=0; i<3; i++) v[i] = B[i] - A[i];
}

// Cross product of vectors made by p1p2,p1p3
// =========================================================================================
// IN:
// p1,p2,p3 - pointers for node coordenates
// OUT:
// normal - pointer for normal vector
void cross(const double *p1, const double *p2, const double *p3, double *normal)
{
	double v1[3], v2[3];
	makeVector(p2,p1,v1);
	makeVector(p3,p1,v2);
	return cross(v1,v2,normal);
}

// Cross product of vectors v1 and v2
// =========================================================================================
// IN:
// v1,v2 - pointers for vectors
// OUT:
// normal - pointer for normal vector
void cross(const double *v1, const double *v2, double* normal)
{
	normal[0] = (v1[1]*v2[2] - v1[2]*v2[1]); // x component
	normal[1] = (v1[2]*v2[0] - v1[0]*v2[2]); // y component
	normal[2] = (v1[0]*v2[1] - v1[1]*v2[0]); // z component
}

// This function is an alternative way to find a tetrahedron vertex opposite to a bdry. face.
// Use FMDB function for this porpouse is memory and cpu costly because user must create
// a structure of faces around tetra. Faces on the interior domain are completely useless!
// Performance test was made comparing local function findVertexOpposite and FMDB equivalent
// and a slighty reduce of time was detected using findVertexOpposite.
void findVertexOpposite(pFace face, pVertex &oppositeVertex)
{	
	list<pRegion> tetraList0,tetraList1,tetraList2;	
	list<pRegion>::iterator iter0,iter1,iter2;
	pPList regionsList;
	pRegion tetra;	

	void *tmp;
	int i;

	set<int> entityNodes;
	vector<pVertex> v;
	M_GetVertices(face,v);
	for (i=0;i<3;i++) 
		entityNodes.insert( EN_id(v[i]) );

	tmp = 0;
	regionsList =  V_regions(v[0]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList0.push_back(tetra);

	tmp = 0;
	regionsList =  V_regions(v[1]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList1.push_back(tetra);

	tmp = 0;
	regionsList =  V_regions(v[2]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList2.push_back(tetra);

	vector<pVertex> tetraVertices;
	bool out0 = true;
	bool out1 = true;
	bool out2 = true;

	for (iter0 = tetraList0.begin(); iter0 != tetraList0.end(),out0; iter0++)
	{
		for (iter1 = tetraList1.begin(); iter1 != tetraList1.end(); iter1++)
		{
			if (*iter0 == *iter1)
			{
				for (iter2 = tetraList2.begin(); iter2 != tetraList2.end(); iter2++)
				{
					if (*iter1 == *iter2)
					{
						tetra = *iter0;
						M_GetVertices(tetra,tetraVertices);
						out0=false;
					}
				}
			}
		}
	}
	tetraList0.clear();
	tetraList1.clear();
	tetraList2.clear();	
	for (i=0;i<4;i++)
	{
		if ( entityNodes.find(EN_id(tetraVertices[i])) == entityNodes.end() )
		{
			oppositeVertex = tetraVertices[i];
			break;
		}
	}
	tetraVertices.clear();
	v.clear();
}

void DijVector(pFace face, pVertex &oppositeVertex, double *Dij)
{
	int i;
	double xyz[3][3], xyzOV[3];

	vector<pVertex> v;
	M_GetVertices(face,v);

	for (i=0 ;i<3; i++) V_coord( v[i], xyz[i] );	
	V_coord( oppositeVertex, xyzOV );

	double a[3] = { xyz[1][0]-xyz[0][0], xyz[1][1]-xyz[0][1], xyz[1][2]-xyz[0][2] } ;
	double b[3] = { xyz[2][0]-xyz[0][0], xyz[2][1]-xyz[0][1], xyz[2][2]-xyz[0][2] } ;
	double c[3] = { xyzOV[0]-xyz[0][0], xyzOV[1]-xyz[0][1], xyzOV[2]-xyz[0][2] } ;

	computeCrossProduct(a,b,Dij);	

	if ( computeDotProduct( c, Dij ) > .0 )
		for (i=0;i<3;i++) Dij[i] = -Dij[i];

	for (i=0;i<3;i++) Dij[i] /= 6.0;
	v.clear();
}

void getTetraWithFace(vector<pEntity> &tetras, pEntity face)
{
	list<pRegion> tetraList0,tetraList1,tetraList2;	
	list<pRegion>::iterator iter0,iter1,iter2;
	pPList regionsList;
	pRegion tetra;	

	void *tmp;
	int i;

	set<int> entityNodes;
	vector<pVertex> v;
	M_GetVertices(face,v);
	for (i=0;i<3;i++) entityNodes.insert( EN_id(v[i]) );

	tmp = 0;
	regionsList =  V_regions(v[0]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList0.push_back(tetra);

	tmp = 0;
	regionsList =  V_regions(v[1]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList1.push_back(tetra);

	tmp = 0;
	regionsList =  V_regions(v[2]);
	while(tetra = (pRegion)PList_next(regionsList, &tmp)) tetraList2.push_back(tetra);

	vector<pVertex> tetraVertices;
	bool out0 = true;
	bool out1 = true;
	bool out2 = true;

	for (iter0 = tetraList0.begin(); iter0 != tetraList0.end(); iter0++)
	{
		for (iter1 = tetraList1.begin(); iter1 != tetraList1.end(); iter1++)
		{
			if (*iter0 == *iter1)
			{
				for (iter2 = tetraList2.begin(); iter2 != tetraList2.end(); iter2++)
				{
					if (*iter1 == *iter2)
					{
						tetra = *iter0;
						tetras.push_back(tetra);
//						M_GetVertices(tetra,tetraVertices);
//						//out0=false;
//						printf("face: %d %d %d\t\t",EN_id(v[0]),EN_id(v[1]),EN_id(v[2]));
//						printf("tetra: %d %d %d %d\n",
//								EN_id(tetraVertices[0]),EN_id(tetraVertices[1]),
//								EN_id(tetraVertices[2]),EN_id(tetraVertices[3]));
//						tetraVertices.clear();
					}
				}
			}
		}
	}
	tetraList0.clear();
	tetraList1.clear();
	tetraList2.clear();
	//exit(1);
}

void getVertexOppositeToFace(pEntity tetra, pEntity face, pVertex &oppositeVertex)
{
	int i;
	set<int> entityNodes;
	vector<pVertex> v;
	
	M_GetVertices(face,v);
	for (i=0;i<3;i++) entityNodes.insert( EN_id(v[i]) );
	
	vector<pVertex> tetraVertices;
	M_GetVertices(tetra,tetraVertices);

	for (i=0;i<4;i++)
	{
		if ( entityNodes.find(EN_id(tetraVertices[i])) == entityNodes.end() )
		{
			oppositeVertex = tetraVertices[i];
			break;
		}
	}	
	entityNodes.clear();
	tetraVertices.clear();
}
